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ABSTRACT 

Molecular clouds are observed to be turbulent, but the origin of this turbulence is not well under- 
stood. As a result, there are two different approaches to simulating molecular clouds, one in which 
the turbulence is allowed to decay after it is initialized, and one in which it is driven. We use the 
adaptive mesh refinement (AMR) code, Orion, to perform high-resolution simulations of molecular 
cloud cores and protostars in environments with both driven and decaying turbulence. We include 
self-gravity use a barotropic equation of state, and represent regions exceeding the maximum grid 
resolution with sink particles. We analyze the properties of bound cores such as size, shape, linewidth, 
and rotational energy, and we find reasonable agreement with observation. At high resolution, the 
different rates of core accretion in the two cases have a significant effect on protostellar system devel- 
opment. Clumps forming in a decaying turbulence environment produce high-multiplicity protostellar 
systems with Toomre-Q unstable disks that exhibit characteristics of the competitive accretion model 
for star formation. In contrast, cores forming in the context of continuously driven turbulence and 
virial equilibrium form smaller protostellar systems with fewer low-mass members. Our simulations 
of driven and decaying turbulence show some statistically significant differences, particularly in the 
production of brown dwarfs and core rotation, but the uncertainties are large enough that we are not 
able to conclude whether observations favor one or the other. 

Subject headings: ISM: clouds - stars:formation - methods: numerical - hydrodynamics - turbulence 



1. INTRODUCTION 

Contemporary star formation occurs exclusively in 
dense molecular clouds (MCs) . Such regions exhibit large 
non-thermal line widths generally attributed to super- 
sonic turbulence (Larson 1981). Although debate contin- 
ues on the origin and characteristics of this turbulence, 
it is now recognized that turbulence is a necessary el- 
ement of star formation and plays an important role in 
the shape of the core initial mass function (IMF), the life- 
times of molecular clouds, and the star formation rate. 

Simulations have shown that supersonic turbulence de- 
cays with an e-folding time of approximately one cloud 
crossing time if there is no energy input to sustain it 
(Stone et al. 1998; Elmegreen & Scalo 2004; Mac Low & 
Klessen 2004). If turbulence decays as quickly in molec- 
ular clouds then star formation must happen rapidly as 
the cloud looses turbulent pressure support and under- 
goes global collapse. In this senario, star formation oc- 
curs on a dynamical timescale and MCs must be transient 
dynamic structures (Elmegreen 2000; Hartmann 2001). 
If however, turbulence is fed from large scales or pro- 
tostellar winds, expanding HII regions, and other pro- 
cesses provide sufficient energy injection to balance dis- 
sipation produced by shocks, then MCs may arrive at 
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a quasi-equilibrium state (e.g., Tan, Krumholz & Mc- 
Kee 2006). Although there are many possible sources 
for turbulent energy the dominant source and the spe- 
cific characteristics of turbulence remain poorly under- 
stood. Recent effort has been directed at this issue and 
some observations of low-mass star forming regions, e.g. 
L1551, find evidence for ongoing turbulence injection 
in the form of winds and jets, which maintains rough 
virial balance over several cloud dynamical times (Swift 
& Welch 2007). While turbulent support is maintained, 
only a small number of overdense regions will become 
gravitationally unstable and form stars in a dynamical 
time, leading to a low star formation rate and allowing 
MCs to live for several dynamical times. 

The two different views of cloud dynamics are related 
to, but distinct from, the two major approaches to sim- 
ulating turbulent molecular clouds. The fact that there 
are two competing approaches to the simulation of such 
clouds is a direct reflection of our lack of understanding 
of the origin of the turbulence in these clouds (McKee & 
Ostriker 2007). In one method, the turbulence is initial- 
ized and then allowed to decay (e.g. Klessen et al. 1998; 
Bonnell et al. 2003; Bate et al. 2002; Tillcy & Pudritz 
2004). The primary problem with this approach is that 
the turbulence decays to levels that are much lower than 
those observed. Advocates of this method argue that 
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the gravitational collapse that ensues after the decay of 
the turbulence can be observationally confused with tur- 
bulence (e.g. Vazquez-Semadeni et al. 2006), but it is 
difficult to see how to maintain a low star formation effi- 
ciency if much of the gas is in a state of gravitational col- 
lapse. The other approach to cloud simulation is to drive 
the turbulence so that it does not decay (e.g. Padoan & 
Nordlund 1999; Gammie et al. 2003; Li et al. 2004; 
Jappsen et al. 2005). This approach allows one to study 
the processes that occur at a given level of turbulence, 
which can be set to match that in a given cloud, but 
it suffers from the disadvantage that the driving is un- 
physical. This approach is a good match for the case of 
quasi-cquilibrium clouds, and it can be made consistent 
with the case of transient clouds if it is assumed that the 
simulation box represents only a small part of the molec- 
ular cloud, so that the decay time for the turbulence is 
long compared to the dynamical time of the simulation. 

The near universal shape of the stellar IMF across di- 
verse star forming environments has sparked much de- 
bate and generated diverse theories. Padoan & Nordlund 

(2004) suggest that the functional form of the IMF can be 
derived from the power spectrum and probability density 
function characteristic of supersonic turbulence. Larson 

(2005) proposes that the peak of the IMF is set by the 
Bonnor-Ebert mass at the minimum cloud temperature, 
which is related to the dust-gas coupling and gas cooling 
efficiency In the competitive accretion model, Bonncll et 
al. (2004) invoke high stellar densities at the centers of 
clusters to propose that the relative position of the stars 
in the gas reservoir determines their mass. According 
to this model, the IMF is determined by mass segrega- 
tion, such that low-mass stars form in the lower density 
gas at the edges of the cluster and higher mass stars 
form in the centers, where their masses can be further 
increased by the coalescence of smaller protostars. In 
addressing the origin of the IMF, numerical simulations 
have been largely inconclusive in discriminating between 
models given that a wide range of conditions (virial pa- 
rameters, resolution, code algorithms, included physics) 
have all succeeded in reproducing the IMF shape. 

A large amount of computational effort has been di- 
rected towards studying self-gravitating turbulent clouds 
both with and without magnetic fields (e.g. Klessen 
2001; Bonnell et al. 2003; Li et al. 2004). A number of 
simulations succeed in reproducing various observed core 
properties such as the IMF and Larson's laws (Padoan 
& Nordlund 1999; Gammie et al. 2003; Tilley & Pudritz 
2004; Li et al. 2004; Jappsen et al. 2005; Bate & Bonnell 
2005). However, most simulations lack the resolution to 
span the turbulent inertial range (Klein et al. 2007) to 
accurately render the evolution of cores into stars in a 
cluster environment. 

In this paper, we perform numerical AMR simulations 
with our code Orion to investigate the role of driven and 
decaying turbulence on low-mass star formation. We fol- 
low the evolution of star forming cores in a turbulent 
box to show that turbulent feedback is correlated with 
the multiplicity of stellar systems, the shape of the IMF, 
and the dominant protostcllar accretion model. In §2, 
we discuss the methodology of Orion and the initial con- 
ditions. In §3, we analyze core properties in driven and 
decaying turbulence at low resolution. In §4, we present 
results from a few high resolution studies of the proto- 



stcllar core evolution inside selected cores formed in the 
context of driven and decaying turbulence. Finally, §5 
contains conclusions. In a companion paper (Offner et 
al. 2008, hereafter Paper II) we investigate the effects of 
driven and undriven turbulence on the properties of the 
cores from which the stars form. 

2. CALCULATIONS 

2.1. Numerical Methods 

Our simulations are performed using the parallel adap- 
tive mesh refinement (AMR) code, Orion, which uses a 
conservative second order Godunov scheme to solve the 
equations of compressible gas dynamics (Truelove et al. 
1998; Klein 1999). Orion solves the Poisson equation us- 
ing multi-level elliptic solvers with multi-grid iteration. 
Throughout our calculations, we use the Truelove crite- 
rion to determine the addition of new AMR grids (Tru- 
elove et al. 1997), 
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where Ax 1 is the cell size on level, I, and we adopt a 
Jeans number of J = 0.25. We insert sink particles in 
regions of the flow that have exceeded this density on the 
maximum level (Krumholz et al. 2004). Sink particles 
serve as numerical markers of collapsing regions and also, 
after sufficient mass accretion and lifetime, protostcllar 
objects. We impose a merger criterion that combines sink 
particles that approach within two grid cells of one an- 
other but prohibits nearby objects from merging if both 
have masses exceeding 0.1 M Q . This limit divides stars 
from substellar objects such as brown dwarfs and has the 
effect of tracking all significantly massive objects. Parti- 
cles that represent temporary violations of the Jeans con- 
dition and have little bound mass tend to accrete little 
and ultimately merge with their more substantial neigh- 
bors. The combination of sink particles and AMR with 
the Jeans criterion allows us to accurately and efficiently 
continue our calculation to high resolution without the 
time constraints imposed by a large base grid size and 
without the consequences of artificial fragmentation. 

2.2. Initial Conditions and Simulation Parameters 

Isothermal self-gravitating gas is scale free, so we give 
the key cloud properties as a function of fiducial val- 
ues for the mean number density of hydrogen nuclei, nu, 
and gas temperature, T. We chose a characteristic 3-D 
turbulent Mach number, M=8A, such that the cloud is 
approximately virialized: 

ha 2 

It is then easy to scale the simulation results to the as- 
trophysical region of interest. For the adopted values of 
the virial parameter and Mach number, the box length, 
mass, and 1-D velocity dispersion are given by 



£ = 2.9 7i 1/2 n Hi 3 /2 pc, 
M = 865 ri 3/2 n Hi 3 /2 Mq , 
ctid = 0.9 Ti 1/2 km s- 1 , 



t ff = 1.37 n^l 2 Myr, 



(3) 
(4) 
(5) 
(6) 
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where n H , 3 = n H /(10 3 cm- 3 ) and T x = T/(10 K) and 
where we have also listed the free-fall time for the gas in 
the box for completeness. For Sh,3 ~ T\ <~ 1, the simu- 
lation approximately satisfies the observed linewidth-size 
relation (Solomon et al. 1987; Heyer & Brunt 2004). For 
the remainder of this paper, all results will be given as- 
suming the fiducial scaling values «h = 1100 cm -3 and 
T=10 K, which are appropriate for the Perseus Molecu- 
lar Cloud (Paper II), but they may be adjusted to differ- 
ent conditions using equations (2)- (5). In terms of the 
Bonnor-Ebert mass, 
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the simulation has a mass of 184M B e- If the Jeans 
mass is defined as Mj = pLj, where Lj = (irc^/Gp) 1 / 2 
is the Jeans length, then Mj = (7r 3 / 2 /1.18)M BE = 

18.97f 2 /<3 M e- 

Our turbulent periodic box study is comprised of two 
stages. The first stage simulates the large scale isother- 
mal environment of a turbulent molecular cloud with 
self-gravity In this "low resolution" stage, we only add 
enough AMR levels to resolve the shape and structure 
of collapsing clumps and cores. This first stage has two 
parts. First, to obtain the initial turbulent spectrum, 
we turn off self-gravity and use the method described 
in MacLow (1999), in which velocity perturbations are 
applied to an initially constant density field. These per- 
turbations correspond to a Gaussian random field with 
flat power spectrum in the range 3 < k < 4 where k is 
the wavenumber normalized to k p i Lys L/2Tr. At the end of 
two cloud crossing times, the turbulence follows a Burg- 
ers P(k) cx k~ 2 power spectrum as expected for hydro- 
dynamic systems of supersonic shocks. For the second 
part, we turn on gravity and follow the subsequent grav- 
itational collapse for two scenarios. It should be noted 
that some workers (e.g. Bate et al. 2003) do not allow 
self-consistent turbulent density fluctuations to build up 
before turning on gravity. Any choice of initialization 
for a turbulent, self- gravitating cloud is necessarily ap- 
proximate, but in our view it is preferable to have self- 
consistent density and velocity fluctuations in the ini- 
tial conditions. In the simulation that we will refer to 
with the letter D (driven), we continue turbulent driving 
to maintain virial equilibrium, while in the other, noted 
with U (undriven), we halt the energy injection and allow 
the turbulence to decay. 

In the second stage, we select a few emerging cores 
for further study in each turbulent box, and we follow 
their fragmentation and evolution into protostellar sys- 
tems at high resolution using a barotropic equation of 
state (EOS). We add additional grid refinement to the re- 
gions we select, which continue to evolve within the low 
resolution context of the box. This method capitalizes 
on the AMR methodology to achieve a high resolution 
study of the development and properties of protostellar 
cores with realistic initial and boundary conditions. Fol- 
lowing all the cores over a free-fall time with AMR rather 
than a subset to the maximum resolution would require 
more than a million CPU hours on 1.5 GHz processors. 
In contrast, our stage two approach with AMR requires 
^50,000 CPU hours per high resolution box. 



In the first stage, it is reasonable to assume that the 
low density gas in the cluster is isothermal and scale- 
free, reflecting the efficient radiative cooling of the gas. 
However, as the gas compresses and becomes optically 
thicker there is a critical density at which the radiation 
is trapped. Ideally, we would directly solve for the ra- 
diation transfer using an appropriate opacity model to 
accurately determine the gas temperature at these high 
densities. However, even approximations such as the Ed- 
dington and diffusion approximations do not sufficiently 
economize the equations of radiation transfer such that 
they are affordable over the resolution and timescales 
necessary for this calculation. Consequently, we adopt 
a bartotropic equation of state to emulate the effect of 
radiation transfer. The gas pressure is given by 



P = P< 



PcC 2 , 
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where c s = (fceT / p) 1 ^ 2 is the sound speed, 7 = 5/3, the 
average molecular weight p = 2.33toh, and the stiffening 
density p c is given by p c /po = 2.8 x 10 s . The value of p 
reflects an assumed gas composition of nuc = O.lnn- The 
value of the stiffening density determines the transition 
from isothermal to adiabatic regimes. It introduces a 
characteristic scale into the previously scale-free isother- 
mal conditions. The isothermal scaling relations above 
remain valid as long as the ratio of stiffening density to 
the average box density is presumed to be constant. We 
chose a value of the stiffening density, p c = 2 x 10~ 13 g 
cm~ 3 , to agree with the p(T) relation calculated by Ma- 
sunaga et al. (1999), who perform a full angle-dependent 
radiation hydrodynamic simulation of a spherically sym- 
metric collapsing cloud core. Unfortunately, we sacrifice 
some accuracy in using the barotropic approximation in 
lieu of radiative transfer, since an EOS assumes that gas 
temperature is a single valued function of density. Simu- 
lations have shown that gas temperature in calculations 
using radiation transfer vs. a barotropic EOS can differ 
by a factor of several and potentially produce different 
fragmentation (Boss et al. 2000; Whitehouse & Bate 
2006). 

The low resolution initial stage uses a 128 3 base grid 
with 4 levels of factors of 2 in grid refinement, giving 
an effective resolution of 2048 3 . The high resolution core 
study has 9 levels of refinement for an effective resolution 
of 65, 536 3 such that the smallest cell size corresponds to 
~10 AU for our fiducial values. 

3. BOUND CLUMP PROPERTIES IN THE 
LOW-RESOLUTION TURBULENT BOX 

3.1. Clump Definition 

At the end of a free-fall time, tg — (37r/32Gpo) 1 ^ 2 , 
with gravity we analyze the core properties and compare 
the driven and decaying turbulent results. At this time, 
the large scale driven turbulent simulation has 32 sink 
particles with 14.2% of the total mass accreted. The de- 
caying turbulence simulation has 20 sink particles con- 
taining 13.6% of the mass. Because the sink particles 
mark collapsing cores rather than individual protostars 
at this stage, these percentages should be viewed as an 
upper limit to the actual star formation rate. Nonethe- 
less, these numbers are not too much larger than the the 
prediction of a 7% star formation rate per free fall time 
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given by Krumholz & McKee (2005) for our assumed 
conditions and neglecting outflows. In the undriven sim- 
ulation, the turbulence decays significantly in ltg and 
no new sink particles are formed after ~ 0.75fg. With- 
out continued driving, there is insufficient energy to cre- 
ate the large scale compressions responsible for seeding 
new cores. After significant turbulent support is lost, 
the cloud deviates from virial equilibrium and the gas 
falls onto existing over-densities rather than forming new 
cores. 

In presenting the results from the low resolution sim- 
ulations, we restrict ourselves to the analysis of ob- 
jects that can best be described as "star-forming bound 
clumps" (see McKee & Ostriker 2007), which are gener- 
ally gravitationally bound but may form several systems 
of stars. In the following sections, we will adopt the ter- 
minology "core" to refer to the bound condensations out 
of which a single protostar (i.e. sink particle) or small 
multiplicity protostellar system forms. Hence, we do not 
apply a Clumpfind algorithm, as described by Williams 
et al. (1994), which also captures unbound and tran- 
sient over-densities. Instead, we define a bound core as a 
sink particle with envelope satisfying four criteria. First, 
the density in the included cells must exceed the average 
shock compressions, i.e., p > p ■ Mm 2 , which also en- 
sure a single peak for each core. Second, the total mass 
in the core must be greater than the local Bonner-Ebert 
mass, signifying that the core will collapse. Each cell, 
i, forming a core must also be individually gravitation- 
ally bound to it such that \E\^\ < |-Ef> E |. Finally, the 
cells included must lie inside a virial radius, R, such that 
avir < 2, where the 1-D velocity dispersion, a, is given by 
the sum of the turbulent and thermal components of the 
gas velocity: a 2 = <j^ t + c 2 . We vary the density cutoff 
by a factor of two and find that changes in the data fits 
remain within 1-sigma error. Thus, our results are not 
overly sensitive to our core definition. The larger of these 
cores may eventually form a cluster of stars and may best 
be described as star-forming clumps. The smaller cores 
will likely form only a single protostellar system. At the 
low resolution stage of analysis it is difficult to predict 
the outcome, and so the line between higher mass star- 
forming clumps and lower mass cores is ill-defined. 

Note that in our methodology, the presence of a sink 
particle does not guarantee the eventual formation of a 
protostar, only that the Jeans condition has been ex- 
ceeded at some time during the simulation. In each sim- 
ulation, there are a few sink particles that do not posses 
envelopes satisfying these criteria. However, all cores in- 
cluded in our analysis are defined to be gravitationally 
bound, collapsing objects rather than transient overden- 
sities in the flow and hence are predisposed to develop 
protostellar systems. 

3.2. Clump Properties 

There are a number of core physical properties that 
are comparable to observations, and we investigate these 
here at 1 iff. Figures 1-5 show the bound core data 
plotted with best fit lines. We exclude objects from the 
fit that have R — y/ab < AAx, where a and b are the 
lengths of the major and minor axes. 
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Fig. 1. — The figure shows the log of the core masses as a function 
of log size (R = \/ab) for the driven (left) and decaying (right) 
boxes at 1 t s . The slopes have fits of 1.03±0.26 and 1.27±0.19, 
respectively. 
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Fig. 2. — The figure shows the core aspect ratios for the driven 
(left) and decaying (right) boxes at 1 iff. The median aspect ra- 
tios for each case are (b/a, c/a) = (0.76, 0.40) and (b/a, c/a) = 
(0.73,0.54), respectively. 
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3.2.1. Density Profiles 
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Fig. 3. — The figure shows the log core velocity dispersions as 
a function of log size (R = \fab) for the driven (left) and decay- 
ing (right) boxes at 1 tg. The slopes have fits of 0.54±0.25 and 
0.19±0.11, respectively. 
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Fig. 4. — The figure shows the rotational parameter, j3, as a 
function of size (R = \/ab) for the driven (left) and decaying (right) 
boxes at 1 tfj. The crosses give the 2-D projected value, while the 
diamonds give the 3-D value. For run D, the median (3 values are 
0.05 (crosses) and 0.05 (diamonds). For run U, the median (3 values 
are 0.08 (crosses), 0.19 (diamonds). 

As plotted in Figure 1, we find that compared to cores 
in run U, the cores in run D have a slightly flatter trend of 
M(R) oc R, consistent with Bonnor-Ebert spheres, which 
are characterized by p(r) oc r -20 . In run U, the cores 
have profiles that are closer to a free-fall profile, where 
p(r) oc r^ 1,5 . Cores that are supported or collapsing 
slowly will tend to resemble pressure-confined isother- 
mal spheres (Kirk et al 2005, Di Francesco et al. 2007) 
as in run D, where turbulence is providing more exter- 
nal pressure support. In run U, where the turbulence 
has decreased significantly, cores tend quickly to infall 
and collapse as unbound gas becomes gravitationally at- 
tracted to the largest over densities. However, the slopes 
of the cores in the two simulations are within 1-sigma er- 
ror due to significant scatter, so that the trends are not 
significantly different. 

3.2.2. Shape 

As shown in Figure 2, both distributions of bound cores 
have similar morphologies and tend to be mainly tri- 
axial. It is thought that in the presence of magnetic 
fields, which we do not include, cores will flatten along 
the field lines (Basu & Ciolek 2004). However, ideal 
MHD simulations by Li et al. (2004) also find that their 
cores are mostly prolate and triaxial 1 . In any event, the 
difficulty of deprojecting observed cores makes the true 
shape distribution ambiguous. Run D has median ma- 
jor and minor aspect ratios of b/a =0.76 and c/a=0.40, 
while the decaying cores have median aspect ratios of 
6/a=0.73 and c/a=0.54. The net medians of the shape 
distributions 0.58 (D) and 0.52 (U) are similar to those 
observed for different star-forming regions which fall in 
the range 0.50-0.67 (Jijina et al. 1999). 

3.2.3. Velocity Dispersion 

In Figure 3, we plot the velocity dispersion as a func- 
tion of core size for comparison against Larson's (1981) 

1 Li et al. 2004 and other references use the word "core" to refer 
to their bound over-densities. For consistency, we continue to use 
our own definition of cores and cores (see § 3.1). 



linewidth-size relation. For low-mass star forming re- 
gions, o~nt R 05 with some sensitivity to core sizes and 
clustering (Jijina et al. 1999). We find exponents of 0.54 
(D) and 0.19 (U). The slope of run D is within the range 
of observed slopes for low-mass regions. Although the 
scatter in our data appears large, our \ 2 fit slope error 
is comparable to the range of fit errors (±0.1 — ±0.19) 
that Jijina et al. report. Plume et al (2000) observed 
massive cores with a completely flat slope, and indeed, 
the cores in run U are more massive with a mean mass 
of 12 M Q versus 8 M Q for the driven, but not significantly 
so (see §3.2.6). A Kolmogorov-Smirnov (KS) test of the 
distributions of velocity dispersions indicates definitively 
that the populations are quite dissimilar at the >> 99% 
level. The difference in slope between the two simula- 
tions is possibly due to crowding in the decaying tur- 
bulent case caused by insufficient global turbulent sup- 
port against gravitational attraction. Jijina et al. (1999) 
showed that clustered objects have a significantly flatter 
linewidth-size relation slope. Note, that the magnitudes 
of the velocity dispersions in run U, although flatter, are 
higher, which is consistent with quickly collapsing rather 
than turbulcntly supported cores. 

3.2.4. Rotation 

Typically, rotational energy makes up only a small 
fraction of the core gravitational energy. The rotational 
parameter (3 is defined as the ratio of the rotational ki- 
netic energy to the gravitational potential energy. For a 
uniform density sphere this can be written: 
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Observationally, fl pos = ft 2 . + fl y is the angular ve- 
locity projected in the plane of the sky, such that 
Prot,obs = |/3 ro t- Goodman et al. (1993), studying 
a selection of dense cores in NH 3 , find that /3 ro t,obs is 
roughly constant as a function of size and find values of 
2 x 10~ 3 < /3 ro t,obs < 1-4 with median /3 rot ~ 0.02. Ob- 
servations of dense cores using N2H + , which primarily 
traces n > 10 5 cm~ 3 , quiescent gas gives similar median 
value of /3rot,obs ~ 0.01 (Casein et al. 2002). For the pur- 
pose of comparison, we evaluate /3 rot in two ways. First, 
we follow the convention of the observers and evaluate 
A-ot.obs by assuming that the cores are projected con- 
stant density spheres. Second we sum over all the 3-D 
data to calculate E rot /E grav . For a singular isothermal 
sphere, E Iot /E grav = |/3 rot . 

Figure 4 confirms that /3 rot for both runs is indepen- 
dent of the core size, and there is fairly large scatter. 
The total range of /3 rot values for observation and sim- 
ulation is roughly the same. We find a range of 0.0005 
< Prot,obs < 0.2 for the driven case and 0.006 < /3 ro t.obs < 
0.3 for the decaying case. However, overall our values are 
a factor of 2 to 4 higher than those found by Goodman 
et al. Run D has a lower median /3 ro t.obs ~ 0.05, while 
run U has very few low /3 ro t,obs cores and so has a me- 
dian /3 ro t,obs ~ 0.08. When we use the complete gas 
properties to calculate E rot /E grav , we find median val- 
ues of 0.05 and 0.19 for the D and U cores, respectively. 
Jappsen & Klessen (2005) perform gravoturbulent driven 
simulations of cores and find a median E Iot /E grav ~ 0.05, 
in agreement with our result. The higher /3 r ot,obs values 
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Fig. 5. — The figure shows the log of the core specific angular 
momentum as a function of log size (R = \/ab) for the driven 
(left) and decaying (right) boxes at 1 iff. The crosses give the 2D 
projected value, while the diamonds give the 3D value. For run 
D, the slopes have fits of 1.91±0.65 (diamond), 1.14±0.31 (cross). 
For run U, the slopes have fits of 1.14±0.35 (diamond), 1.50±0.23 
(cross) . 

measured in the cores in the undriven simulation may be 
a side effect of the smaller turbulent support: Since the 
U cores are moving more slowly, they may more easily 
accrete gas from farther away, which has higher angular 
momentum. (We thank the referee for this comment.) 
Although a KS test verifies that the two /3 ro t,obs popula- 
tions are distinct, neither is a good match for observation 
since both have median values that arc higher than ob- 
served. 

One possible explanation for the factor of 3-5 differ- 
ence between simulation and observation is that mag- 
netic fields play a significant role in decreasing core rota- 
tion. A number of recent simulations of isolated rotating 
magnetized cloud cores have shown that magnetic brak- 
ing is an efficient means of outward angular momentum 
transport (Hosking & Whitworth 2004; Machida et al. 
2004; Machida et al. 2006; Bannerjee & Pudriz 2006). 
The oblate cores formed in the ideal MHD simulation of 
Li et al. (2004) show a median /3 ro t,3D similar to ours 
(Li, private communication), however, all their cores are 
supercritical by an order of magnitude. 

Another possibility to account for the difference in me- 
dian (3 is that observers typically investigate isolated 
cores, which are easier to distinguish and analyze but 
tend to be less turbulent. However, our study specifically 
concerns bound cores forming in a turbulent cluster. Us- 
ing Larson's laws, /3 rot oc v 2 ot /(GM/R) oc R/(GM/R) oc 
1/S ~const. However, there is large scatter and a few 
exceptions of clouds with non-constant column density, 
S, such that measurements of (3 ro t could be sensitive to 
differences in column density in various MCs. 

3.2.5. Angular Momentum 

There is also a substantial difference between the spe- 
cific angular momentum in the two cases as illustrated 
in Figure 5. We plot both the 3D total specific angular 
momentum of the cores, which is obtained by directly 
summing the angular momentum of the individual cells 
comprising a clump, and the 2D specific angular momen- 
tum, by totaling the projected momentum along a line 



of sight. In run D, the specific angular momentum fits, 
J2d(-R) oc R and J3d(-R) oc R , bracket the expected 
j(R) oc R 15 based upon the linewidth Sv oc R 1 / 2 and 
assumption of virial balance (Goodman et al. 1993; this 
argument suggests the same value for both the 2D and 
3D cases). The specific angular momentum fits in run 
U are more similar but still a little flat, j2T>{R) oc R 15 
and J3d(-R) oc R 11 . Because the decaying cores are less 
turbulent, they will be inclined to have less variation of 
angular momentum than their driven counterparts. The 
cores in run D overshoot the expected relationship for 
j 3D , while the decaying cores undershoot by a similar 
amount. In either case, we expect that the simulated 
angular is affected by the absence of braking effects from 
magnetic fields, which we do not include in these sim- 
ulations. Nonetheless, we find that the measured range 
of j <~ 10 21 — 10 22 cm 2 s _1 to be consistent with ob- 
servational estimates and the 2D angular momentum es- 
timates to be statistically similar to one another, but 
flatter than the measured j2D oc R 1 - 6 ^ - 2 (Goodman et 
al. 1993; see Fig. 5). 

3.2.6. Core Mass Function 

Measurements of the core mass function (CMF) show 
that its shape strongly resembles the stellar initial mass 
function (IMF) (Lada et al. 2006). The high mass end, 
in particular, seems to share a similar power law index. 
As a key characteristic of star formation, the core and 
star mass functions for driven and undriven turbulence 
has been extensively numerically studied. Ballesteros- 
Paredes et al. (2006) and Padoan et al. (2007) find a 
mass function of the form dM /dlog(m) oc to -1 ' 3 for cores 
in driven hydrodynamic turbulence, even without the 
presence of self-gravity. Klessen (2001) finds that both 
driven turbulence with 1 < k < 2 and undriven turbu- 
lence produce a core spectrum with a similar slope to that 
of the measured IMF. A number of isothermal SPH sim- 
ulations of decaying turbulence have shown agreement 
with the observed IMF despite different initial turbu- 
lent conditions in which a turbulent velocity spectrum is 
initialized on a constant density field and then allowed 
to decay in the presence of self-gravity (e.g. Klessen & 
Burkert 2001; Bate el al. 2002; Bonnell et al. 2003; Tillcy 
& Pudritz 2004; Bonnell et al. 2006). In this method, 
the turbulence does not reach a steady state and the 
simulated cloud is not virialized as observed. 

For the purpose of comparison, we plot the CMF for 
the simulations D and U at 1 tg in Figure 6. The two 
runs produce 30 and 19 bound cores, respectively. Un- 
fortunately, the statistics at the high mass end are too 
small to be able to rule out either distribution on the ba- 
sis of agreement with the Salpeter slope. Although the 
agreement looks better for the driven cores, we find that 
the mass distributions are in fact statistically similar ac- 
cording to a KS test. 

Overall, the simulations have statistically different dis- 
tributions of angular momentum, rotational parameter, 
and velocity dispersion. The decline of turbulent com- 
pressions in the undriven run appears to make some sig- 
nificant changes and causing fewer new condensations to 
be formed. As turbulent pressure support is lost, the 
contracting gas instead falls onto existing cores resulting 
in less turbulent, more quickly rotating cores than in the 
driven case. However, it is not possible at this time to 
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Fig. 6. — The figure shows the sink (dashed line) and core (solid 
line) mass distributions for the driven (left) and decaying (right) 
runs at 1 tg. The straight time has a slope of -1.3. 

say which approach corresponds more closely with obser- 
vation. 

4. PROTOSTELLAR CORES AT HIGH RESOLUTION 
4.1. Overview 

In this section, we present our computational results 
for the evolution of the protostellar systems contained 
in a few selected cores using 5 or 6 additional levels of 
refinement. We accomplish our study by inserting a re- 
finement box around the core of interest before a sink 
particle is introduced on level 4 so that cells inside the 
box continue to higher densities and refine according to 
the Jeans criterion, while cells in the remainder of the 
simulation refine only to a maximum level of 4 as be- 
fore. The lengths of the high resolution boxes are typ- 
ically 0.25-0.5 pc depending on the size of the enclosed 
clump. As a result, the boxes contain a region ^200-2000 
times volumetrically smaller than the simulation domain. 
The total initial mass in the boxes ranges from 4-12 M Q , 
which easily encompasses the bound core and all collaps- 
ing regions associated with it. In this way, we can per- 
form a high resolution study of selected collapsing cores 
with realistic initial conditions and consistent boundary 
conditions taken from the surrounding lower resolution 
grids computationally cheaply and efficiently. In each 
portion of the box, sink particles are introduced when 
the corresponding maximum refinement level is reached. 
At high resolution, each sink particle represents a single 
"protostellar core." Thus, we are able to follow the clump 
fragmentation at high resolution, without the need to re- 
run the entire calculation at that resolution. 

We chose six cores for further study. In cases Ula 
and Ulb, we test for convergence by following the same 
cores at two different resolutions. In cases D2 and U2, 
we choose an early collapsing object that is present in 
both the driven and decaying simulations to highlight 
the differences between the calculations. The cores Ul, 
D2, and U2 have initial bound masses greater than 2.5 
M . We also study two smaller driven and undriven 
cores, D3 and U3. As a result, we observe the effect of 
turbulent support on protostellar system development. 
The physical properties of the selected cores are given in 
Table 1 . Table 2 gives the three dimcnsionless quantities 
relating to the box surrounding each core: M, a vir , and 



the self-gravity parameter, 



// : 



M 



where a v ; r can then be written as 



M 2 



2/3 



(10) 



(11) 



These three parameters characterize the amount of 
turbulence in the core vicinity, the degree of self- 
gravitization of the gas, and the extent to which balance 
is achieved between the two. Table 2 indicates that all 
the small boxes are subsonic and thus the influence of 
gravity is dominating the gas in the regions around the 
cores. 

Note that when we cease driving in run U, the turbu- 
lent cascade continues and the turbulent decay rate is 
determined by the Mach number and the domain size as 
described by Mac Low (1999). At any given time, the 
effect of the decay on the cores forming in the high reso- 
lution subdomain depends upon the amount of turbulent 
decay in the large box. The 1-D velocity dispersion in 
Table 1 is an indicator of the change in turbulent energy 
when the core of interest is collapsing. 

4.2. Convergence Study 

Before embarking on further analysis, it is important 
to show that the results at the calculation resolution are 
suitably converged. In particular, it is necessary to show 
not only that there is no artificial fragmentation but that 
the number of fragments is constant with increasing res- 
olution. For our convergence study we consider a box in 
the decaying turbulence run, Ul, which encloses a long 
filament that collapses to form a number of small over- 
dense fragments along its length. We run this calculation 
with 9 (Ulb) and 10 (Ula) levels of refinement, which 
corresponds to a minimum cell size of ^10 and 5 AU, 
respectively. Figures 7, 8 and 9 show the two simula- 
tions at 16 kyr, 23 kyr, and 53 kyr, respectively, after 
the formation of the first sink particle. Tables 3 and 4 
give the sink particle masses and the fragment masses at 
these times. We define the fragments as discrete cores of 
bound gas with density greater than 2 x 10~ 16 g cm~ 3 . 

We find that both resolutions produce the same num- 
ber of collapsing fragments and yield a similar collection 
of sink particles. Most of the fragment masses for the 
two resolutions differ by at most a few percent, while 
sink particle masses may differ by 50%. At a particular 
instant in time, discrepancies between the number of sink 
particles in the two runs can occur due to several factors. 
First, the addition of extra levels allows the higher res- 
olution simulation to collapse for a longer time without 
exceeding the Jeans criterion. Thus, a sink particle ul- 
timately forms in both cases at the same location but 
at slightly different times. Another possibility is that a 
sink particle forms in both simulations at similar loca- 
tions, but in one it mergers with a larger neighbor. A 
final possibility is that the region that collapses in the 
higher resolution becomes thermally supported before a 
sink is formed. In all these cases the gas physics can 
be quite similar but the introduction of sink particle can 
differ due to small details. For example, at 3 kyr the low 
resolution simulation has formed sink particles in each 
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TABLE 1 

LOW RESOLUTION CORE PROPERTIES FOR EACH CASE WHEN THE SINK PARTICLE HAS 0.1 Mq . 





Ula/Ulb 


D2 


U2 


D3 


U3 


Core Mass (Mq) 


10.71 


5.05 


4.32 


1.59 


2.80 


imax (PC) 


0.45 


0.23 


0.19 


0.03 


0.08 


Shape 


1:0.28:0.05 


1: 0.37:0.09 


1:0.66:0.26 


1:0.69:0.60 


1:0.78:0.24 


U rms (km/s) 


0.42 


0.32 


0.36 


0.33 


0.45 


n avc (10 5 cm- 3 ) 


1.44 


1.08 


1.01 


4.06 


1.25 


Ovir 


1.5 


1.12 


1.90 


1.06 


2.24 


/3rot 


0.052 


0.025 


0.011 


0.028 


0.047 


tff(ioV) 

M 1U a 


8.9 


10.3 


10.7 


5.3 


5.9 


2.4 


1.9 


3.9 


4.9 


3.0 



a m lD i s the velocity dispersion of the entire box, which is fixed at 4.9 for the driven cases. 



TABLE 2 

Turbulent box properties for the whole domain and the 
small boxes containing the cores. 





D/U at t=0 


Ula/Ulb 


D2 


U2 


D3 


U3 




8.37 


0.71 


0.80 


0.72 


0.38 


0.74 


Mbox 


206.82 


19.03 


10.74 


12.41 


2.6 


4.59 


Ct vir 


1.67 


0.06 


0.11 


0.08 


0.06 


0.16 



Note. — The values for the small boxes are determined using 
the length L sma n=0.25 pc. 



filamentary fragment with condensation masses ranging 
from 1.5 x 1CT 2 - 8 x 1CT 2 M Q (see Table 4), while the 
higher resolution run has not reached sufficient density 
for any sink particles to form. This rather odd filamen- 
tary structure is created and confined by the ram pres- 
sure of intersecting shocks. As a result it forms somewhat 
smaller bound clouds than the minimum Bonnor-Ebert 
mass associated with the local pressure (P~3x 10 6 dy 
cm~ 2 ) at p ~ 1CP 14 g cm~ 3 . The smallest sink particles 
formed in the filament later merge as shown in Figure 
8 when the gas in the filament streams onto the disk- 
protostar system. 

At later times and for small masses the corresponding 
sink particle properties differ the most significantly, par- 
ticularly at the lower mass end as shown in Table 3. Due 
to the intrinsically chaotic and dynamically unstable na- 
ture of three or more body systems, at later times the 
evolution of the two calculations begins to diverge. This 
is unsurprising because not only do the calculations have 
different AMR grid structures, but the particle members 
of the system are introduced at slightly different times 
and initial masses. Despite this, the masses and config- 
uration still show reasonable agreement at 80 kyr. 

4.3. Influence of Turbulence on Stellar Properties 

Interstellar turbulence undoubtedly has a substantial 
effect on cloud lifetimes and core creation, however, its 
relationship with core fragmentation and evolution is less 
certain. The level of turbulence in cores is partially de- 
pendent on how much mass and energy the envelope ex- 
changes with the surrounding turbulent gas. In turn, 
the properties of the parent core influence the rate of 
protostellar core formation and accretion. If substantial 
mass continues to fall onto the clump, as in the case of 
global contraction, then external flow patterns will im- 
pinge upon on the system development, increasing the 
accretion rate and possibly causing fragmentation. If 



TABLE 3 

Masses of the stars in Mq for decaying simulations at two 
different resolutions at two different times after the 
formation of the first sink particle. 



At 23 kyr 53 kyr 



Resolution 5 AU 


10 AU 


5 AU 


10 AU 


0.834 


0.705 


1.224 


0.925 


0.000 


0.216 


0.000 


0.369 


0.264 


0.262 


0.571 


0.455 


0.171 


0.175 


0.762 


0.768 


0.000 


0.036 


0.106 


0.141 






0.061 


0.036 






0.128 


0.180 



Note. — The subscripts 10 (Ula) and 9 (Ulb) represent the 
number of AMR levels. The sink particle absence in the second 
row of the high resolution column is due to an early merger (m j 
O.IMq with the neighbor listed in the first row. 
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Fig. 7. — The figure shows the log column density of a core in 
the decaying turbulence simulation Ulb (left) and Ula (right) with 
resolution of 10 AU and 5 AU 16 kyr after the formation of the 
first sink particle. 

however, the core accretes at a relatively low level in the 
manner of Bondi-Hoyle accretion in a turbulent medium 
then the core will accrete much less over time (Krumholz 
et al. 2006; Krumholz et al. 2005). Protostars forming 
in such a core limit to the Bondi-Hoyle accretion rate as 
the high density gas is depleted. 

In cases D2 and D3 we continue turbulent driving to 
maintain virial equilibrium. 

To avoid directly adding artificial perturbations that 
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TABLE 5 

Masses of the protostars in both driven and decaying 
simulations at 260 kyr (the larger cores, d2 and u2) and 
130 kyr (the smaller cores, d3 and u3). 



540 

X[AU] 



Fig. 8. — The figure shows the log column density of a core in 
the decaying turbulence simulation Ulb (left) and Ula (right) 23 
kyr after the formation of the first sink particle. 



-0.33 0.33 




-16201 

-1620 -810 810 1620 
X[AU] 



10 810 1620 
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Fig. 9. — The figure shows the log column density of a core in 
the decaying turbulence simulation Ulb (left) and Ula (right) 53 
kyr after the formation of the first sink particle. 



TABLE 4 

Core gas mass (Mq) including embedded sinks for the 
decaying simulations at two different resolutions at 
three different times. 



A t 


16 kyr 


23 kyr 


53 kyr 


Resolution 


5 AU 


10 AU 


5 AU 


10 AU 


5 AU 


10 AU 




0.741 


0.758 


0.997 


0.997 


1.907 


1.905 




0.124 


0.124 


0.280 


0.283 


0.827 


0.806 




0.077 


0.077 


0.177 


0.180 


0.136 


0.191 




0.040 


0.037 












0.035 


0.033 












0.291 


0.259 




















0.106 


0.036 



Note. — The subscripts 10 (Ula) and 9 (Ulb) represent the 
number of AMR levels. The minimum density of the gas is p = 
2 X 10 — 16 g cm -3 . The represent cores that have merged with 
others and cannot be individually distinguished. 



D2 (Mq) 


U2 (M Q ) 


D3 (Mq) 


C/3(M ) 


1.221 


1.811 


0.639 


0.586 


1.047 


1.002 


0.453 


x 0.552 


1.049 


0.933 




0.348 


0.490 


x 0.223 




x 0.114 


x 0.382 


x 0.131 




0.048 


0.329 


0.059 




x 0.047 


0.281 


x 0.034 






0.207 


x 0.030 








x 0.023 




a 



Note. — The x's represent particles that are ejected from the 
system by dynamical interactions. The time of first sink particle 
formation after the onset of gravity for each of the cores is 270, 
680, 250, and 660 kyr for D2, D3, U2, and U3, respectively. 

"One additional BD mass sink particles have been excluded as 
numerical disk fragmentation from column U2. 

may affect the core development or seed new fragmenta- 
tion, we do not apply any velocity perturbations to the 
high resolution regions inside the refinement box. Thus, 
the turbulence cascades into the highly-refined box from 
the outside in a self-consistent manner. In cases U2 and 
U3, the simulation continues without any turbulent in- 
jection. 

We find striking differences in the protostellar systems 
formed in the driven and decaying cores. The most obvi- 
ous difference between the D and U runs is the difference 
in the number and mass of sink particles formed (Table 
5). For example, initially the fragmentation of D2 and 
U2 is similar temporally and spatially, but U2 eventually 
forms a slightly larger number of objects particularly at 
small masses as the level of turbulence in the two sim- 
ulations diverges. The small D3 core forms a small sta- 
ble binary system at early times, whereas U3, which is 
also fairly small, fragments into a number of protostellar 
members. 

Despite the small-number statistics, we are able to 
compare the IMF of the protostars to the observed initial 
mass function via the KS test. The KS test determines 
the probability that a given data set is drawn from a 
specified statistical distribution, in this case, the single 
star IMF given by Chabrier (2005) . This test is accurate 
for input sets of 4 or more data points. We achieve a 
best fit by scaling the masses by an adjustable normal- 
ization, e = m^/m s - m ^, where m s i n ]< is the mass of the 
sink particle. Given that this simulation lacks feedback 
effects such as outflows and radiation transfer, the sink 
particle masses represent an upper limit and the scaling 
factor corresponds to an efficiency factor of e=0.25-0.75 
(Matzner & McKee 2000). 

Figure 10 shows the scaled cumulative distribution 
function (cdf) for the runs D2, U2, and U3, and the cdf 
of the Chabrier is overlaid for comparison. Although all 
three runs can have a high confidence level of agreement 
with the measured IMF, the normalization values and 
the shapes of the distributions are quite different. For 
example, the smaller stellar population of D2 has fewer 
low-mass objects and hence has a smaller efficiency scal- 
ing factor of ~ 0.4 with highest likelihood of being drawn 
from the IMF of 67%. Conversely, U2+U3 distribution 
contains collections of low-mass objects and intermediate 
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Fig. 10. — The figure shows the cumulative distribution function 
(solid line) at t=0.26Myr for D2 (left), U2+U3 (rig/ii),where the 
dotted line is the Chabrier 2005 IMF fit. The dashed vertical line 
represents the point of largest disagreement. The probability that 
the data are drawn from the Chabrier IMF is 67% and 59%, re- 
spectively, where the efficiency scale factors of the simulations are 
0.4 and 1.0, respectively. 

mass objects, where the largest disagreement occurs in 
the middle of the two populations. A scaling factor of 1.0 
gives the best agreement of 59%. A scaling factor near 
unity implies that protostellar mass loss has a negligi- 
ble effect on the final mass of the star, contrary to some 
theoretical expectations (Shu ct al 1987; Nakano et al. 
1995; Matzner & McKee 2000). For the D2 distribution, 
the largest disagreement occurs at the higher mass end, 
indicating that if the protostars continue to accrete mass 
and no new protostars are formed, then it is likely that 
the high probability of agreement with the IMF will be 
maintained while the scale factor shifts to a lower value. 
The U2+U3 have a larger scale factor due to the signif- 
icant number of low mass objects with accretion halted 
by dynamical ejection. These objects will be unlikely 
to accrete additional mass and are essentially fixed. For 
the undriven runs, the largest disagreement occurs in the 
middle of the distribution, indicating a widening differ- 
ence between the sub-stellar fixed-mass ejected objects 
and those that remain in the gas reservoir and continue 
accreting. Further running time will more likely make 
the gap wider and agreement worse. 

The efficiency scale factor is also dependent upon the 
normalization we have chosen. The minimum mass that 
we are able to resolve in these simulations is proportional 
to the Jeans mass evaluated at the maximum level of 
refinement. Because the Jeans mass is inversely related 
to the density, normalizing the results to a density higher 
than our fiducial value of n = 1100 cm~ 3 will produce 
lower mass objects and shift the IMF peak towards lower 
mass. This will also increase the efficiency factor used to 
scale the distribution to the universal IMF. 

Studying the time evolution of the two simulations 
shows the origin of the different stellar populations. In 
simulation D2, the initial collapse and core fragmentation 
produces three well separated objects that remain fairly 
far apart. A few additional objects form, but they do 
not suffer large gravitational interactions with the pri- 
maries and so remain in the high density regions and 
continue accreting. By contrast, in U2 and U3 the lack 



of global turbulent support causes mass to fall onto the 
early formed protostellar cores resulting in contraction 
of the clump. This causes all the protostars to gravitate 
towards the core center. As the protostellar proximity 
increases, the accretion disks interact and become grav- 
itationally unstable (see discussion in § 4.6). Fragmen- 
tation ensues. The stellar systems become increasingly 
dynamically unstable with the addition of these small 
latecomers, which rapidly suffer strong gravitational in- 
teractions with the larger protostars and are thrown out 
of the high density reservoir of gas. Their small envelopes 
are stripped away, thus truncating the accretion process 
and effectively fixing their stellar masses (see Figure 11). 
This truncation process occurs for approximately half of 
the objects formed in the undriven simulations (see dis- 
cussion of brown dwarfs in § 4.5) 

4.4. Accretion 

There are two main accretion paradigms. In both mod- 
els, star formation begins with the outside-in formation 
of gravitationally bound cores and their inside-out col- 
lapse. However, the core accretion model proposes that 
the main protostellar accretion phase takes place early 
on and continues until the entire core mass is accreted or 
expelled (Shu et al 1987; Nakano et al. 1995; Matzner 
& McKee 2000), after which accretion becomes negligi- 
ble. Thus, the initial core size and subsequent feedback 
effects limit the mass of the protostars. In contrast, the 
competitive accretion model proposes that stars begin in 
a core as wandering, accreting 0.1 M Q seeds, whose fi- 
nal mass is determined by the protostar's location in the 
clump (Bonnell et al. 1997, 2001). Mass segregation is 
a common feature of this model, such that the largest 
mass objects inhabit the region of highest gravitational 
potential and the smallest objects inhabit the less dense 
gas, usually having been ejected from the center by grav- 
itational interactions. 

In our results, there are two accretion phases. Initially, 
there is a transient period of high accretion during which 
the initial infalling gas accretes onto the newly formed 
sink particle and gas is depleted from the cells inside 
the accretion region. This phase is model independent 
and occurs while the newly created sink particle region 
reaches pressure equilibrium with the surrounding gas. 
Generally less than 10% of the accretion occurs during 
this time. During the second phase, the accretion rate 
approximates the Shu model for core accretion, 

M. = c 3 JG, (12) 

(Shu et al. 1987), although in most cases the accretion 
rate is gradually declining. This solution is valid until 
the rarefaction wave reaches the core, i.e. when approx- 
imately half of the original core mass has been accreted 
(McLaughlin & Pudritz 1997; McKee & Tan 2002). After 
that, the accretion rate diminishes as the density of the 
surrounding gas decreases, but our simulations end be- 
fore it is possible to determine if this final stage actually 
occurs. 

To illustrate the differences between the protostellar 
systems, we have plotted the mass as a function of time 
for all sink particles, the instantaneous mass accretion 
rate for the first two formed objects, the time-averaged 
accretion rate, and the total mass in sink particles as 
a function of time. The turbulent core accretion and 
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Fig. 11. — The figures show the sink particle mass as a function 
of time for runs D2, D3, U3, and U2 shown clockwise from top left. 
Each particle is represented by a different style line. 
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Fig. 12. — The figures show the instantaneous sink particle 
accretion rate as a function of time for runs D2, D3, U3, and U2 
shown clockwise from top left. Only the history of the two first 
forming particles is shown. 

competitive accretion models describe the evolution of 
the stellar population in cases D2 and U2, respectively. 
In the former case, objects are mainly formed from core 
fragmentation with separations larger than 1000 AU. The 
average accretion onto the protostars initially agrees with 
the Shu model but, modulo fluctuations, diminishes over 
time as the core mass depletes (Figure 12). Meanwhile, 
the core envelope accretes according to the Bondi-Hoyle 
model of turbulent accretion (Krumholz et al. 2006). For 
driven turbulent environments the overall Mach number 
will remain sufficiently high such that the core will not 
gain a substantial amount of mass during the core dy- 
namical time and the main accretion phase of the form- 
ing protostars will be limited by this time. However, 
in the decaying turbulent case loss of turbulent pressure 
support potentially causes significant additional mass to 
accrete onto the core, resulting in a more constant proto- 
stellar accretion rate (Figure 12, bottom row). However, 
the differences in the accretion rates of the most massive 
objects are subtle due to the significant fluctuations. 

Perturbations to the accretion disks and dumpiness of 
the infalling gas cause fairly large variability in the sink 



Fig. 13. — The figures show the averaged sink particle accretion 
rate for the first two sink particles as a function of time for runs 
D2, D3, U3, and U2 shown clockwise from top left. The average is 
taken over 10 consecutive timesteps, and the solid flat line indicates 
the value of cf/G. 

particle accretion rate as illustrated in Figure 12. How- 
ever, we do not observe that most of the mass is deposited 
in short intervals by dumpiness in the disk as noted by 
Basu et al. (2006), who model 2D axi-symmetric disks 
with magnetic fields. The absence of this effect in our 
calculations is most likely due to our Cartesian geometry 
rather than lack of magnetic fields (Basu, private com- 
munication). The r — 4> geometry used by Basu et al. is 
more suitable for disk treatment and has lower numerical 
viscosity, which may suppress small scale dumpiness. 

Due to differences in core accretion, the two cases pro- 
duce much different stellar populations. In the driven 
cases, in which fewer objects form, protostars accrete 
more smoothly and do not undergo strong dynamical in- 
teractions with their neighbors. However, in the decaying 
cases, the strong infall pushes the protostars to the core 
center and the large number of nearby objects accret- 
ing from the central reservoir of gas causes the smallest 
objects to be kicked out of the cluster. This can be ob- 
served in Figure 12 (U3) as a precipitous drop off in the 
accretion rate for individual objects or as flatlining of 
the object mass (Figure 13, bottom row). Differences be- 
tween the number of objects are caused by turbulent sup- 
port, which prohibits new mass from infalling onto the 
clump. High accretion causes fragmentation and drives 
forming objects to the gravitational center, where they 
dynamically interact. This close proximity results in ob- 
ject ejection and destabilization of the accretion disks, 
which leads to new fragmentation. The differences in 
accretion rate and stellar population between the two 
cases suggests that the maintenance of turbulence, pro- 
tostcllar accretion, and stellar population arc intimately 
related (Krumholz et al. 2005). In spite of individual 
accretion fluctuations, the total accretion of the objects 
from the core is dominated by the largest objects, such 
that the fraction of the core accreted is relatively smooth 
over time as shown by Figure 14. 

4.5. Brown Dwarfs 

Brown dwarfs (BD) are observed to comprise ~ 10-30 % 
of all luminous objects in star forming regions (Andersen 
et al. 2006; Luhman et al. 2007). For example, in the 
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Fig. 14. — The figures show the total mass accreted normalized 
to the initial bound core mass as a function of time for runs D2, 
D3, U3, and U2 shown clockwise from top left. 

Chabricr (2005) IMF, with which we compare, BDs with 
masses M* < 0.08M Q comprise ^20% of the total num- 
ber of objects. Understanding the population, origins, 
and connection between planets and hydrogen-burning 
stars is essential in formulating a successful theory of 
star formation. Observations remain particularly am- 
biguous concerning the primary formation mechanism of 
BDs, sparking many theories. Of these, proposals for BD 
formation by turbulent fragmentation, ejection, or via 
disk fragmentation have the most potential for generat- 
ing BDs in sufficient numbers (Padoan & Nordlund 2004; 
Reipurth & Clarke 2001; Whitworth et al. 2007). Sim- 
ulations provide an important vehicle for testing these 
theories, and we discuss the BD population of our simu- 
lations in this section. 

Our driven turbulence simulations D2 and D3 do not 
produce any sink particles of final substcllar mass, which 
are primarily created in our simulations by prematurely 
truncated accretion. However, if BDs are formed via 
the same mechanism as stars, accrete from a disk, and 
produce outflows (Luhman et al. 2007) then the same 
efficiency factor will be valid for scaling all the objects in 
the simulation. The decaying turbulence simulations U2 
and U3, however, produce a much larger initial number 
of BDs, 33%. This agrees with the competitive accretion 
paradigm: Bate et al. (2002) find that ~44 % of the 
objects that form in their 50M Q decaying turbulent cloud 
simulation qualify as BDs, which is much higher than the 
observed fraction. Their calculation is initialized with 
uniform density and an initial turbulent velocity field, 
however they do not drive the turbulence, so that the 
turbulence never achieves a steady relaxed state. Despite 
this difference, their result is in qualitative agreement 
with the BD formation mechanism and number fraction 
of our decaying turbulence runs. 

Ideally, we would like to understand the BD popula- 
tion in various star forming regions as a function of their 
general properties. The turbulent fragmentation model 
predicts an upper limit on the total mass available for 
the formation of BDs as a function of the Mach number 
and average number density (Padoan & Nordlund 2004) . 
According to their model, the total gas mass available 
to make BDs from turbulent compressions is 0.4% of the 
total gas mass or 3.7 M as a function of our simulation 



Mach number and density. If the SFR per free fall time 
for the driven and decaying runs is respectively 14.3% 
and 13.6% then the total maximum possible mass in BD 
due to turbulent fragmentation as a fraction of the actual 
mass turned into stars is 2.8% and 2.9%. For compari- 
son, the fraction of the actual luminous mass turned into 
BDs according to the Chabrier IMF is ~2%. Our high 
resolution protostellar systems have a BD mass fraction 
of 0.0% and 3.2% for D2, and 1X2 + U3, respectively, 
using the efficiency factor from Figure 10 and includ- 
ing M* < 0.08-Mq in the BD mass total. However, the 
turbulent fragmentation model gives only the maximum 
fraction of gas that can be converted to BDs by turbulent 
compressions and it does not include possible BD forma- 
tion in disks (Goodwin & Whitworth, 2007) . Fragmenta- 
tion of disks and dynamical ejection is responsible for all 
of the BDs in the decaying simulation. Thus, comparison 
between the decaying turbulence models and turbulent 
fragmentation prediction is misleading. 

The absence of BDs in the driven runs is reasonable if 
BDs actually form via turbulent fragmentation. In such 
a process, small low-mass objects form from small low- 
mass cores. Since we have not chosen any particularly 
small cores for high resolution study, we would not ex- 
pect to find many BDs. Thus, scaling to the stellar IMF, 
which has a peak at <~ 0.2 M Q requires a small efficiency 
factor. The core distribution is set by the low resolution 
turbulent initial conditions. The minimum expected core 
mass is the Bonnor-Ebert mass evaluated at the maxi- 
mum turbulent gas density. According to the turbulent 
fragmentation model, this maximum is set by the prob- 
ability density function (PDF) of the gas density. The 
resolution and Mach number of our simulation yield a 
density PDF that falls off at - 1.3 x 10~ 18 g cm" 3 or 
Mbe — O.2M . The minimum mass estimated from this 
density agrees with the minimum sink particle mass at 
the end of a free fall time at low resolution. Since this 
mass is well above the maximum BD mass, it also ex- 
plains the low abundance of low-mass objects at high 
resolution in the driven simulation. Moreover, this sug- 
gests that the driven high-resolution IMF distribution is 
incomplete at the low-mass end such that scaling to the 
actual IMF may be optimistic and result in underesti- 
mating the core efficiency factor. 

One observational measure of the BDs in a region 
is given by the ratio of low-mas stars to BDs: R = 
N(0M - 1.0M Q )/(AT(0.02 - 0.08). Measurements of lo- 
cal star-forming regions give a range of R ~ 2 — 5 (e.g. 
Andersen et al. 2006). For the driven and decaying 
simulations, respectively, we find R > 7 and R = 3.0, 
although these numbers are clearly sensitive to the low- 
statistics of our simulations. These ratios are most prop- 
erly represent lower limits because we have not included 
the effects of radiative transfer, which have been shown to 
suppress fragmentation (Krumholz et al. 2007). Over- 
all, the driven high-resolution BD mass fraction more 
closely agrees with the turbulent fragmentation predic- 
tion, whereas the undriven BD mass fraction agrees more 
closely with competitive accretion results, which predict 
larger BD numbers and mass fractions. Given that BD 
formation via disk fragmentation dominates in the un- 
driven case, it is unsurprising that these statistics do not 
agree well with the turbulent fragmentation model and 
are quite different from one another. 
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4.6. Disk Stability 
Analytically, gravitational disk instability is dictated 



by the Toomre Q parameter, which is given by 

Q= "" 



ttGS : 



(13) 



where n is the epicyclic frequency, and £ is the surface 
density. For values of Q < 1, the disk becomes unsta- 
ble to gravitational fragmentation. Spiral arms develop 
for low Q values and fragmentation ensues when Q ap- 
proaches 1 from above. This fragmentation manifests as 
a density increase at those locations. The early frag- 
mentation in D2 and U2 generally occurs near the disk 
perimeters, where it is coldest (e.g. Figure 15). In the 
simulations, sources of disk instability are due to a combi- 
nation of perturbations from clumpy infalling gas, grav- 
itational influence of nearby bodies (i.e. other sink par- 
ticles), and from actual collisions between disks. Since 
the sink accretion radius is r acc = AAx, we neglect the 
innermost 4 cells in the analysis. We define the disk gas 
where p > 2 x 1CP 16 g cm~ 3 , which agrees fairly well 
with disk boundaries determined visually. In general, we 
find disk radii between 150-300AU and surface densities 
of a few g cm~ 2 , values similar to observed properties of 
low-mass disks (Andrews & Williams 2006). 

The stability of a disk and the onset of gravitational 
instability have been shown to be correlated with the 
accretion rate of the disk itself (e.g. Bonnell 1994; Whit- 
worth et al. 1995; Hennebelle et al. 2004; Matzner & 
Levin 2005). Higher disk accretion rates increase the 
likelihood of disk instability. This fact agrees with our 
observation that more instances of disk fragmentation 
occur in simulations U2 and U3, where there is larger in- 
fall onto the clump, in contrast to the case D2 where the 
disks remain fairly stable. The level of disk instability is 
directly visible in the plots of the sink particle accretion 
rates (Figure 12); very noisy and irregular accretion cor- 
responds to clumping and disturbance of the disk. The 
simulations where sinks have many close neighbors show 
the highest rates of disk instability and episodic accre- 
tion. Note that in Figure 12 of run U4, the ejection of a 
companion substantially reduces the accretion rate fluc- 
tuations of the remaining protostar. 

There has been considerable recent discourse on the 
necessary criteria for resolving disks and preventing ar- 
tificial fragmentation (Nelson 2006; Klein et al. 2007; 
Vorobyov & Basu 2005; Durisen et al. 2007). Since we 
do find that our disks fragment, this is a topic of con- 
cern. Most recent simulations, including ours, have used 
the Jeans condition as defined by Truelove et al. (1997) 
or Bate & Burkert (1997) to set the minimum refine- 
ment of meshes and particles, respectively, in the disk 
under study. However, Nelson (2006) argues that this 
criterion is inadequate and inappropriate for cylindrical 
disk geometry. Additional possible sources of error in 
our calculation may arise from sink particle gravitational 
softening, numerical viscosity, and the cartesian nature 
of the AMR grid. We address these issues here. 

In calculating the gravitational sink particle-particle 
and sink particle-gas interactions, we use a constant soft- 
ening length 0.5Ax max , where Ax max is the grid spacing 
on the maximum level. This is much smaller than both 
the disk radius and the size of the accretion region, so 
it should have little effect on the behavior of the disk. 
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Fig. 15. — The figure shows the log column density (g cm -2 ) 
of an accretion disk in run D2 with 11 levels of refinement. Two 
fragments have formed at the edges of the spiral disk structure. 
The solid black lines denote grid boundaries, a 

In general, we find that the observed disk fragmentation 
occurs as cores forming in the ends of spiral arms, well 
away from the center of the disk (see Figure 15). 

Nelson (2006) requires two specific criteria for adequate 
disk resolution. The first is a Toomre condition, 



T> 



Ax 1 
At 



where T is the Toomre number, 
stable wavelength defined by: 



(14) 

and At is the neutral 
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and Ax 1 is the cell spacing on level, I. The above criterion 
is analogous to the Jeans criterion defined in Truelove et 
al., 

Ax 1 

~a7' 



j > 



(16) 



For our simulations, a disk radius of 200 AU is cov- 
ered by 20 or 40 cells, which is fairly marginal resolu- 
tion, but we will show it is, in fact, sufficient. We plot 
the azimuthally averaged density and Toomre Q param- 
eter (equation 13) as a function of radius in Figure 16. 
Density enhancements are correlated with low Toomre 
Q in each refinement case. We also plot the right hand 
sides of equations (14) and (16) as functions of radius in 
Figure 17. In all cases, these quantities are under the 
fiducial value of 1/4. The over resolution in the central 
part is due to the requirement that all cells surrounding 
a sink particle be refined to the maximum level in order 
to encompass the accretion region. 

Figure 15 indicates the borders between AMR grids, 
so some disk regions within 200 AU become derefined, 
and Ax 1 — > Ax 1 ^ 1 . However, these regions still satisfy 
the refinement criteria by a good margin. 

The second criterion formulated by Nelson applies 
specifically to SPH codes. It postulates the necessity 
of resolving the disk scale height at the midplane by four 
smoothing lengths. Nelson argues that insufficient res- 
olution of the vertical structure produces errors in the 
force balance, thus favoring artificial fragmentation. If 
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we assume a one-to-one conversion between smoothing 
lengths and grid cells, we can apply it to our calcula- 
tion. Figure 17 shows azimuthally averaged quantities 
for an accretion disk for ~2.5, 5, and 10 AU maximum 
resolution. The lowest resolution run fails to adequately 
resolve the disk scale height, but we do not see extra 
fragmentation. We may not see this effect because Nel- 
son formulated and tested his criteria for SPH rather 
than grid-based codes. It is also possible that the one- 
to-one analog of smoothing length to Ax is not the cor- 
rect conversion. However, most disagreement is in the 
inner regions where the artificial viscosity is high, which 
suppresses any potential fragmentation. In order to de- 
termine the cause of the discrepancy, further exploration 
with a full grid high resolution study of disks is necessary. 

A careful study of the sink particle accretion in a 
disk is given in Krumholz et al. (2004). For a Kep- 
lerian disk, Orion agrees well with analytic results ex- 
cept for some irregularities when the radius of the disk, 
r ~ r B — GM/c^, the Bondi radius. However, our sim- 
ulations have r << re during the main accretion phase 
and should be unaffected. Also of concern is the magni- 
tude of the numerical viscosity, which has the potential 
to suppress fragmentation if it is sufficiently high. Using 
the definition of a viscosity defined in Krumholz et al. 
(2004) , we can estimate the magnitude of this viscosity 
as a function of disk radius: 

^fJLY 385 (17) 



a~78 



Ax \AxJ 
-O.SM^Ax, 
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where is the radial distance from the central star in 
units of 150 AU, Mi is the stellar mass is units of M©, 
Ax$ is the cell size in units of 5 AU, normalized to the 
maximum level of refinement, and T w is the gas tem- 
perature in units of 10K. This expression shows fairly 
large sensitivity to the cell size and disk radius. Due to 
the large a value in the inner region of the disk, artifi- 
cial viscosity is likely to significantly influence the disk 
properties within the inner 100 AU. Large values of a 
could potentially suppress disk fragmentation. Proto- 
stcllar disks around low-mass protostars, which are fairly 
thin and have a low ionization fraction, have been mea- 
sured to have viscosities of a ~ 0.01 (King et. al 2007; 
Andrews & Williams 2006). 

We find that all disks form exactly two fragments at 
the same radial locations where Q ~ 1. Convergence of 
the disk density distribution and number of fragments is 
our main concern. The averaged quantities are slightly 
different in the three cases, although the general trends 
are the same. In the lower resolution case the fragmen- 
tation is less pronounced, however two sink particles are 
introduced at these locations. It is certainly true that the 
fragments are not well resolved at the lowest resolution, 
and we are only marginally resolving the disks. Hence 
we do not devote much discussion in this paper to disk 
properties. Serious study of disks requires much higher 
resolution than we adopt in this paper and is best studied 
in cylindrical or polar coordinate geometry to minimize 
the effects of Cartesian cell imprinting. 

Given that observations find stellar systems have typi- 
cally 2-3 stars (Goodwin & Kroupa 2005), the large num- 
ber of objects produced in the high resolution decaying 
simulations seems somewhat anomalous. However, the 
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Fig. 16. — The figure shows azimuthally averaged disk properties 
as a function of log radius (AU) for a disk with ~ 2.5 AU (top), 
5.0 AU (middle) and 10.0 AU (bottom) resolution. The left plots 
shows log p for the edge on view of the disk. Plots on the right 
show log Q vs. log r. The central region corresponding to the sink 
particle accretion region is excluded from the plot. 

actual initial multiplicity is more difficult to determine 
than multiplicity among field stars due to the difficulty 
of detecting small obscured objects, some of which may 
have separations below the resolvable limit. In addi- 
tion, systems with more than two bodies are unstable 
and decay through gravitational interactions ultimately 
decreasing the multiplicity of stellar systems over time. 
We witness this behavior in the decaying turbulence pro- 
tostellar systems, which expel low-mass members. 

Nonetheless, it is likely that a few of these small frag- 
ments are numerical products, resulting from our equa- 
tion of state. For example, Boss et al. (2000) and 
Krumholz et al. (2006) both find that fragmentation 
is sensitive to thermal assumptions and the inclusion of 
radiative transfer, since heating tends to enhance disk 
stability. Price & Bate (2007) show that magnetic fields 
tend to suppress and delay both fragmentation and spiral 
disk structure. It is probable that inclusion of radiative 
feedback and magnetic fields would suppress some of the 
small objects that we find in the undriven runs. However, 
the absence of these objects in the driven simulations in- 
dicates a striking difference in the accretion rate, system 
stability, and fragmentation history with turbulent feed- 
back. 

5. CONCLUSIONS 

In this paper we use turbulent simulations with AMR 
to illustrate distinctions between driven and decaying 
turbulence. Despite identical initial conditions in the 
two simulations, we find significant differences between 
the two cases after one free-fall time. Our simulations ne- 
glect the effect of magnetic fields, which are poorly obser- 
vationally constrained and occupy a place of ambiguous 
but potentially large importance (Crutcher 1999). Our 
simulations also lack radiation transfer, instead relying 
on the barotropic approximation, which may affect the 
core fragmentation and protostcllar multiplicity in our 
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Fig. 17. — The figure shows azimuthally averaged disk properties 
as a function of log radius (AU) for a disk with ~ 2.5 AU (top), 5.0 
AU (middle) and 10.0 AU (bottom) resolution. The first column 
shows plots of J (dashed line) and T(solid line) vs. log r, where 
the horizontal line marks the fiducial value of 0.25. The second 
column shows the number of cells in the disk vertical scale height 
as a function of log r. The solid line is the required resolution of 
the vertical scale height according to the Nelson criteria and the 
dashed line is our resolution. The central region corresponding to 
the sink particle accretion region is excluded from the plot. 

results. 

We find that the properties of the cores in driven and 
decaying turbulence at low resolution are not sufficiently 
different to completely dismiss one turbulent environ- 
ment. This is in part due to the large scatter in our 
results. For example, we find that the cores in the differ- 
ent environments have similar shapes and mass-size rela- 
tions. However, we find that cores in the driven simula- 
tion have less rotational energy, which is in better agree- 
ment with observations (Goodman ct al. 1993; Caselli et 
al. 2002). The linewidth-size relation of the cores form- 
ing in driven turbulence is also closer to the observed re- 
lation for low-mass regions (e.g. Jijina & Adams 1999), 
while the linewidth-size relation of cores in the decaying 
simulation is quite flat. We find that driven turbulence 
produces a greater number of cores than decaying turbu- 
lence with the potential for new star formation occurring 
longer than a single dynamical time. In contrast, the de- 
caying simulation stops forming new condensations be- 
fore one global free-fall time. 

The largest differences between the two cases are ap- 
parent at high resolution. We show that our high reso- 
lution simulations are converged and that the resolution 
is sufficient to capture core fragmentation, despite being 
marginal for determining the detailed properties of disks. 
We find that the presence or absence of global virial bal- 
ance has only a subtle influence on individual accretion 



rate of the largest object forming in the core at least 
for the first few core free fall times. However, the cores 
forming in a decaying turbulence environment show clear 
signs of competitive accretion such that a core's accretion 
rate is tied to its dynamical history and and its location 
in the clump. This supports the results of Krumholz et 
al. (2005) who show that simulations exhibiting com- 
petitive accretion do so because of lack of a source of 
turbulence. 

The loss of turbulent feedback in the decaying run af- 
fects the dynamic behavior of the forming protostars, 
resulting in significant disk fragmentation, and BD for- 
mation by ejection. This leads to overproduction of BDs 
in comparison to the observed IMF (e.g. Chabrier 2005). 
In contrast, the driven simulations form few BDs, which 
can be understood in the context of the turbulent frag- 
mentation model for star formation, which predicts BDs 
to mainly form from small highly compressed cores. Ob- 
servations of low-mass star forming regions do not find 
large velocity differences or significant spatial segrega- 
tion between BDs and low-mass objects as obtained in 
the decaying simulation. 

While our simulations of driven and decaying turbu- 
lence show some statistically significant differences, par- 
ticularly in the production of brown dwarfs and core ro- 
tation, the uncertainties are large enough that we are not 
able to conclude whether observations favor one or the 
other. However, in Paper II we use simulated line profiles 
to estimate core velocity dispersions and centroid veloci- 
ties, and we find that decaying turbulence leads to highly 
supersonic infall onto protostars, which has not been ob- 
served. Our results thus give some support to the use 
of driven turbulence for modeling regions of star forma- 
tion, but a conclusive determination of which approach 
is better awaits larger simulations with the inclusion of 
magnetic fields, protostcllar outflows and thermal feed- 
back. 
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